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We present a method for computing the evolution of a spacetime containing a massive particle 
and a black hole. The essential idea is that the gravitational field is evolved using full numerical 
relativity, with the particle generating a non-zero source term in the Einstein equations. The matter 
fields are not evolved by hydrodynamic equations. Instead the particle is treated as a quasi-rigid 
body whose center follows a geodesic. The necessary theoretical framework is developed and then 
implemented in a computer code that uses the null-cone, or characteristic, formulation of numerical 
relativity. The performance of the code is illustrated in test runs, including a complete orbit (near 
r = 9M) of a Schwarzschild black hole. 

PACS numbers: 04.25.Dm, 04.20.Ex, 04.30.Db, 95.30.Lz 
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I. INTRODUCTION 

This paper is concerned with the evolution of a space- 
time containing a small object in orbit near a black hole, 
the goal being to compute the motion of the object and 
the emitted gravitational radiation. In the techniques 
normally used for this type of 'radiation-reaction' prob- 
lems, the small object is treated as a point-particle evolv- 
ing on a fixed background spacetime with its self-force 
taken into account. There are a number of approaches 
concerning the implementation of the self- force - for ex- 
ample [JSI!ll,|Ell0SHm[H[l2 A n alternative 
approach that could be used is full numerical relativity 
including computational relativistic hydrodynamics (see 
for example ES El El El El El , 

and for a review 191). 
To our knowledge, such simulations have not been per- 
formed for the extreme mass ratios considered in this 
paper. 

The approach developed here uses full numerical rela- 
tivity but with the hydrodynamic aspect greatly simpli- 
fied, and can be described as a polytropic particle (PP) 
method. We use full numerical relativity for the evolution 
of the gravitational field (including a non-zero stress en- 
ergy tensor as source) . The matter fields are not evolved 
by relativistic hydrodynamics but rather the object is 
treated as a quasi-rigid body whose center is evolved by 
the geodesic equation. This approach avoids the intrica- 
cies and computational expense of relativistic hydrody- 
namics. 

Of course, there are limitations to the PP approach: it 
can only be applied in situations in which the internal hy- 
drodynamics of the material object are unimportant, and 
in which the rigidity approximation is reasonable. For ex- 
ample, one would certainly need full numerical relativity 
and relativistic hydrodynamics in situations in which the 



object is expected to be tidally disrupted. On the other 
hand, if tidal effects are small (in a sense that will be 
made precise), then the PP model should provide a good 
description of the physics. 

The bounds of the domain of applicability of the PP 
method have not been investigated systematically. How- 
ever, from the example runs presented in Sec. the 
method should be applicable under appropriate restric- 
tions, e.g. when m < 10~ 5 M where M is the mass of the 
black hole and m is the mass of the particle, and when 
the size of the particle is small enough that it not be 
tidally disrupted. Thus, we expect the particle method 
to be applicable in astrophysical situations involving the 
inspiral and capture of a neutron star or white dwarf by 
a galactic black hole (with M about 10 6 M Q ). On the 
other hand, it is difficult to see how the method could 
be used for a stellar remnant black hole (with M about 
WMq) - any object with small enough m/M would have 
too large a diameter. Thus the method is expected to 
make predictions concerning gravitational radiation that 
will be relevant to observations by LISA [2(j , rather than 
to observations by LIGO or other earth based detectors. 

The results presented here use a characteristic gravity 
code 0,H2|. However, there is no reason to restrict the 
PP method to the characteristic approach. The method 
should also be applicable within Cauchy formulations of 
numerical relativity. 

We have also constructed a finite difference version of 
a 8— function model. We found that the PP model per- 
forms better, giving smoother results and, in particu- 
lar, exhibiting convergence with grid-refinement (as de- 
scribed in Sec. IV A|l . For these reasons, although we de- 
scribe both models, we give implementation details and 
results only for the PP model. 

The purpose of this paper is to introduce the PP 
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method. We also give some example runs exhibiting in- 
spiral and a plunge to the black hole. These runs demon- 
strate the potential of the method, and are not an ac- 
curate description of the physics - in particular, in the 
inspiral case, the inspiral rate is much larger than that 
predicted by other methods. Of course, validated phys- 
ical results are the goal of this project but, as discussed 
in Sec. lVll more computational testing and development 
is needed before the goal can be attained. 

We begin by summarizing previous results on the char- 
acteristic formulation of numerical relativity in Sec. [H] 
Issues concerning the theoretical framework of a massive 
particle are discussed in Sec. IIIII Sec. IIVI presents, in 
detail, the computational algorithms. Tests of the code 
and example runs are given in Sec. E] 



The Einstein equations decompose into hypersurface 
equations, evolution equations and conservation laws. 
The hypersurface equations form a hierarchical set for 
v >r , fc r , /3 r , B r , (r 2 Q) ir , U >r and W >r ; and the evolu- 
tion equation is an expression for (rJ). ur - The explicit 
form of the equations is given in |22j in the vacuum case. 
The matter source terms, in the case of a perfect fluid of 
density p, pressure P and velocity v a , are stated in [32| : 
except the matter source term in Eq. H2-(31) is incor- 
rect and that equation should read 



2(rJ), ttr - ((l + r- 1 ^) (rJ) r ) ^ 
+2r- 1 e /3 3 2 e /3 - (r -1 W) r J + Nj 
4e 2/ V(p + P) 



H7) 



((JV ang - KV ang f + V 2 ng ) , (2.4) 



II. SUMMARY OF PREVIOUS RESULTS, AND 
NOTATION 



The formalism for the numerical evolution of Einstein's 
equations, in null cone coordinates, is well known pH 
El E3 (see also jM El El El El) For the sake of 
completeness, we give here a summary of the formalism, 
including some of the necessary equations. The version of 
the gravity code being used here is fully described in E3 • 

We use coordinates based upon a family of outgoing 
null hypersurfaces. We let u label these hypersurfaces, 
x A (A = 2, 3), label the null rays and r be a surface area 
coordinate. In the resulting x a = (it, r, x A ) coordinates, 
the metric takes the Bondi-Sachs form 28, 301 



ds = — 



w 

(! + -)" 



2r 2 h AB U B dudx A 



r 2 h AB U A U B j du 2 - 2e 2p dudr 

„2 J. . _ j„Aj„B 



r hAsdx dx 



(2.1) 



where h AB hsc = $c an< ^ det(fiAB) — det(qAB), with 
qAB a unit sphere metric. We work in stereographic co- 
ordinates x A — (q,p) for which the unit sphere metric 
is 



q A Bdx A dx B = -j=^(dq 2 



dp 2 ) , where F = 1 + q 2 



\-p- 

(In previous notation we used P = 1 + q 2 + p 2 . Here we 
change notation to F because P now represents pressure, 
which we cannot denote by p because that is a stereo- 
graphic coordinate). We also introduce a complex dyad 
q A = -j(l,z) with i = \/— 1. For an arbitrary Bondi- 
Sachs metric, Hab can then be represented by its dyad 
component 



J = h AB q A q B /2 



(2.3) 



with the spherically symmetric case characterized by J = 
0. We introduce the (complex differential) eth operators 
3 and S (see [3lJ for full details), as well as a number of 
auxiliary variables K — h,AB<l A q B '/2, U — U A qA, Qa = 



r 2 e- 2 ^h A BU B , Q 



-- h AB q A q B /2, U 
) A q A , B = 3/3, v = 9J and k 



where V ang = VAq |32(, and Nj is defined in |2JJ 
and |22| . The remaining Einstein equations reduce to 
conservation conditions which need only be satisfied on 
the inner boundary, which are automatically satisfied 
here because the boundary has a simple Schwarzschild 
geometry. 

The null cone problem is normally formulated in the 
region of spacetime between a timelike or null worldtube 
r and 1 + . We represent I + on a finite grid by using 
a compactified radial coordinate x = r/(l + r). The 
numerical grid is regular in (x, q,p) and consists of two 
patches (north and south), each containing n x n q n p grid- 
points. The x— grid covers the range [0.5, 1]. Each angu- 
lar grid patch extends two grid-points beyond the domain 
(q,p) G [-q s ,qs] x [-q s ,qs], with q s > 1. Thus there is 
an overlap region at the equator with larger overlap for 
larger q s . 

We denote the Bondi-Sachs metric (|2.1|l by g a p and 
the background metric {g a p with J = U = = 0, W = 
—2M/r) by g\M]a/3- The mass M of the black hole is 
normally scaled to M = 1 in simulations. 



III. THEORETICAL FRAMEWORK 

We have developed two different particle models with 
rather different conceptual frameworks but implemented 
with very similar numerical codes. This section describes 
each of the two frameworks, as well as some other theo- 
retical issues. 



A. The polytropic particle (PP) model 

The PP model treats the particle as an object of fixed 
size in its local proper rest frame, with its center z a (u) 
describing a geodesic of the full spacetime. The particle 
is treated as quasi-rigid, and a simple formula is used to 
evaluate the stress-energy tensor, which then appears as 
source in the Einstein equations. The simplest model of 
a polytrope, which will be used here, is for the case with 
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index n = 1. Then the density p and pressure P of an 
equilibrium configuration are given by |33j| 



B.TT 



P 



157 o 2i?^ : 



2^2 



4.R.R 2 



P 



(3.1) 



for R < i?*, and p = P = for R > where i?» is 
the radius of the polytrope and R is the distance from 
its center. 

The density p(R) and pressure P(r) at a point x a = 
(u, x l ) in the Bondi-Sachs coordinate system are set by 
determining the distance i?, in the local proper rest 
frame, between x a and the geodesic described by the 
center of the polytrope. We first define the displacement 
vector e Q relative to the polytrope's center, 



e u = 0, 



z\u). 



(3.2) 



The projection of e a into the hypersurface orthogonal to 
the worldline at time u is 



R a = {gap + v a vp)e 13 , 



(3.3) 



with g a p evaluated at z Q , and v a evaluated at time u. 
Then we define the proper distance R as the magnitude 
of R a , 



R = Jg^RaRp, 



(3.4) 



and set the density and pressure using Eq. I|3.1|l . Thus, 
in its local proper rest frame, a (small) polytropic particle 
is a spherically symmetric object with proper radius 
The perfect fluid condition is used to set the stress-energy 
tensor, i.e. T a/3 = (p + P)v a Vp + Pg a p. 

The PP model is approximate in the following sense. 
Equation l|3.1[) is exact only for an isolated sphere in 
equilibrium in Newtonian theory. In general relativity 
it is a good approximation if m -C R* (and the example 
runs presented later satisfy this condition). Further, in 
the field of a black hole, a polytrope would not preserve 
its spherical shape but would become tidally distorted. 
Both these approximations introduce the same type of 
error - the PP's motion is treated as quasi-rigid but, 
for the purpose of determining its gravitational field, its 
stress-energy tensor T a p is modeled as a perfect fluid. 
Thus there are contributions to T a p that are being ig- 
nored, and the magnitude of these contributions is now 
estimated. 

Firstly, since the simulations in Sec. satisfy the con- 
dition m <C R* , the error in ignoring the tidal stress will 
dominate that of using Eq. (|3.1I) to estimate the pressure. 
We use Newtonian theory to make an order of magnitude 
estimate of the tidal stress. The tidal acceleration in the 
radial direction is 2Mx/r 3 where x is the distance from 
the center of the polytrope. The tidal stress is maximum 
at the center of the polytrope, and a simple dimensional 
argument (which can be confirmed by integration) shows 
that 



Cmax — O 



Mm 



(3.5) 



In order to estimate the significance of Umax in the stress- 
energy tensor, we compare it to the maximum density 



fmax 
Pmax 



= s = o 



ARiM 



(3.6) 



In a physical situation in which the polytrope is sta- 
ble against tidal disruption, (MRl)/(mr 3 ) < 0(1), and 
the ratio S is smaller than 0(mj i?*), which, as previ- 
ously noted, is small here. However, in the example to 
be considered in Sec. IV Bl the polytrope is not stable 
against tidal disruption, and thus the internal stresses 
may be large. In fact, the ratio S is approximately 
4/(81tt) « 0.016. 



B. The S function model 

At the analytic level, a point particle of mass m at 
position z a = (it, z % ) = (u, r, z A ) has density p and 4- 
velocity v a satisfying 



p^f^gv 11 = m5(x l - z l ), 



(3.7) 



where J 8(x l — z l )drdqdp = 1. We model the (5-function 
on the grid by assigning weights w to each grid point 
in a stencil surrounding the particle. In terms of a test 
function <f>, this requires 



■>i<Pi 



W]A\ 



(3.8) 



where is a sum over a stencil I of grid points sur- 
rounding the particle position z l and Ay is the coordi- 
nate 3-volume of the stencil /. We determine the weights 
Wi representing the (5-function by choosing a set of test 
functions, e.g. for the stencil of eight points determined 
by the cell surrounding the particle we choose 

(p = a + a,i(x l — z l ) + a id (x l — z l )(x 3 - z j ) 

+ aijkix* - z*)^ - z j ){x k - z k ), (3.9) 

where i ^ j ^ k so that the a's constitute eight arbitrary 
coefficients. This then gives 8 simultaneous equations to 
solve for the wj, which are given explicitly in Eq. I|4.6|) 
below. 

It is necessary to renormalize the metric so as to avoid 
infinities in the equations of motion. The metric occurs 
through the normalization of the 4-velocity v a and the 
raising of indices. We take the components (v r , Va) to be 
basic since they represent the pullback of the 4-velocity 
to the null hypersurface. We renormalize the other com- 
ponents by using the background metric g\M]ap to raise 
indices and to normalize the 4-velocity. This avoids the 
problem of an infinite self-potential energy of the particle 
and is in keeping with the principle that the energy of the 
particle only depends on its velocity and position in the 
Schwarzschild field. It should be emphasized that this 
renormalization, or use of g[M] a p rather than g a p, ap- 
plies only to the undifferentiated metric. Metric deriva- 
tives that occur in the particle equations of motion are 
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computed using the full metric g a @ - otherwise radia- 
tion reaction could not be included and we would sim- 
ply be computing the motion of a test particle in the 
Schwarzschild geometry. Of course, it is the full metric 
which is evolved by the characteristic algorithm. 

Although we were able to use the <5-function model to 
compute qualitatively reasonable orbits for the problems 
considered in Sees. IV Bl and IV 01 they were significantly 
less smooth than with the PP model, and the growth 
in the deviation from a Schwarzschild orbit was much 
faster. Further, we did not find any quantity that ex- 
hibited convergence, making it problematic to use the 
5- function model to obtain physical predictions. Perhaps 
a more sophisticated numerical approach might lead to 
convergence of global quantities, but this would be a dif- 
ficult project that we do not pursue here. 



C. Modeling the particle orbit 

A goal of this work is to study radiation reaction. 
This is a small effect, and, numerically, it could be hid- 
den if terms of order unity (representing the background 
Schwarzschild geometry) are added and subtracted in the 
equations of motion. The motion of a test particle in 
a background Schwarzschild geometry satisfies certain 
conservation laws; the motion of a particle with mass 
< m < M does not satisfy these laws, but they are 
nevertheless useful because these laws indicate quanti- 
ties that change very slowly. Our strategy is to use the 
Schwarzschild conservation laws to find quantities that, 
in the general case, evolve slowly. In the process, all 
background Schwarzschild terms cancel out and we are 
left with expressions involving only small quantities. 

For the case of a test particle in the Schwarzschild ge- 
ometry there is a reflection symmetry plane, the plane 
of the orbit. Thus the normalized 4- velocity is com- 
pletely determined by its components T a v a = v u and 
$> a v a = Vcj), where T a and <& Q are the Killing vectors of 
the Schwarzschild background and is a velocity com- 
ponent with respect to (u, r, 9, </>) null-spherical coordi- 
nates. In the general case, v u and are approximately 
conserved. Partly because of the stereographic coordi- 
nates being used, the implementation is quite technical 
(see Sec. USEl for details). 



D. Caustics 

The characteristic evolution code breaks down if caus- 
tics develop, which render the null coordinate system 
employed singular. A rough estimate can be readily ob- 
tained by employing the well-known condition for the 
deflection of light by a massive body such as the Sun. 
We find as an approximate condition for caustics not to 
form that 

-g > r. (3.10) 

4m 



IV. COMPUTATIONAL METHOD 
A. Overview 

The PP method evolves both the matter and gravity 
fields. At each grid-point at which the density is non- 
zero, the particle's density, pressure and 3-velocity are 
found and used to construct the right hand side of the 
Einstein equations, and the gravitational field is then 
evolved as described in Sec. |nj The gravitational field 
affects the motion of the particle: the 3-velocity v% is 
evolved by using the geodesic equation in the form 



dvj _ T mi vsv t g aS g ie % 
du v u 



(4.1) 



(4.2) 



and the particle's position is evolved by 

dz l v l 
du v u 

The setting of initial data is described in Sec. II V CI 
below. The worldtube T at r = 2M is the (past) hori- 
zon of a Schwarzschild black hole of mass M . Thus the 
boundary data on T has the simple analytic form |32| 

J = is = k = f3 = B = U = Q = 0, W = -2M. (4.3) 



B. Computational algorithms 

The iterative evolution algorithm proceeds as follows: 

1. Start at time u = u^ n \ The gravitational field 
5o/3 _1 ' > ^ s k nown over the whole grid and the bound- 

(n) 

ary data supplies g a i in a neighborhood of r = 2M . 

The particle's position z 1 ^ and velocity are 
also known. 

(n) 

2. Determine the grid-cell G P containing the point 
z l (™); i.e., determine a, such that 



„(ai) 



< z 



l(n) 



< r 



(ai + l) 



q M <z 2(n) <g (a 2 + l) 5 
r,(«3) < z 3(n) < (a 3 + l) _ 



(4.4) 



This is done on both north and south patches, al- 
though if the particle is not in the equatorial over- 
lap region there will be a solution for only one of 
the patches. We define 

S Mn) _ z i{n) _ x ai (4 5) 

and we define weights at the eight grid-points at 
the corners of Gp ^ by 

Mn) A 2(n) c3(n) 
/ ai+ei \ _ °ei »e 2 °e 3 

w[x > ~ A 1 (»)A 2 WA 3 W 
where e, = {0, 1}. 



(4.6) 
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3. Next, we set the density and pressure. In general, 
this needs to be done on both north and south 
patches. The density at the grid-point x l is set 
by means of Eqs. I|3.1[l - il3.4|) ; then the pressure P 
is set. 

4. The Einstein equations are now integrated to find 
the metric g^ . The source terms are given in |32T ] 
(where, as already noted, Eq. [§^-(31) should be 
replaced by Eq. (J23J). 



5. The formula g a ^v a vp 



(n) 

- 1 is used to find vi. ; 



the 



metric g a @ is known at the required grid-points and 
its value at the particle position z^ 71 ' is found by 
taking a weighted average using the weights found 
in Eq. ( 14. 6f) above. 

6. The formula v a — g a ^vp is used to find v aty7l \ again 
using the weighted average to find g a @ . 

7. Equation l|4.2[l is now used to find z ^ n+1 \ On the 
first time-step, this is done by the Euler method 
and, on subsequent time steps, by the 3-point 
Adams-Bashforth method, i.e. 



Au / dz 
~2~\ du 



We now find u t - +1 \ The right hand side of 
Eq. (|4.1|l is evaluated at z^ n ' by, as usual, finding 
the value at the grid-points and taking a weighted 
average using the weights found in Eq. f !4.6l) above. 
The terms in Eq. I|4.1|l are quite complicated and 
were found using a Maple script, which was also 
used to generate Fortran code. Details are given in 
an Appendix. The numerical evolution method is 
the same as used in step[7| 



(n) 



dz_ (n - 
du 



(4.7) 



C. Setting the initial data 

We have experimented with various ways of setting the 
initial data, and found that, apart from some early tran- 
sient effects, the options tried make little difference to 
both the particle orbit and the gravitational field. The 
initial gravitational content is prescribed by setting J = 
at u = Q. (It is also possible to set the initial data J by a 
Newtonian limit condition p3ll34ll35| . the computational 
implementation of which will be discussed elsewhere). Of 
course, this means that we are introducing spurious gravi- 
tational radiation into the initial data, which might have 
the effect of initializing the particle into a different or- 
bit to that intended. We tried evolving the code for a 
pre-determined time us (and us could be possibly set to 
zero), during which time the particle's velocity and posi- 
tion are not updated. The idea is that the gravitational 
field should relax to the correct form, with the spurious 



initial gravitation content radiating away, by the time us 
when the particle is allowed to move. 

The code requires the initial velocity as a 1-form Vi 
but a physical description normally specifies the tangent 
vector v % . For example, a particle in a circular orbit 
would have 



v r = 0, (v p ) 2 + (v q f 



F 2 M 



4r 2 (r-3M)' 



(4.8) 



Suppose that we are given v l rather than Vi. Initially, 
when only the background metric is known, Vi is con- 
structed from v* using 



-1 



(4.9) 



to first determine v u ; then Vi = g[M]ia vC *- Then, while 
u < us, the code uses the fact that v a is found at each 
time step to determine a value of v r such that v r = 
by an iterative algorithm. Explicitly, we use the secant 
algorithm 



v r ( a +i) = v r (a) - v r (a) fs 



' ] r{a) 



V 



r(a-l) 



(a) (a-1) 



(4.10) 



where a is the iteration number and fg is a factor (which 
is 1 in the standard algorithm) that may need to be set 
to 0.2 or smaller for stable convergence - the difficulty 
here is that we are solving v r (v r ) = not as a simple 
algebraic equation but as an equation whose coefficients 
change as the metric relaxes. 

We found that the value of us has little effect on com- 
puted orbits, at least for the cases computed in Sec. IV Bl 
below. Thus, we will set us — 0. Nevertheless, the option 
of setting a non-zero Ms is retained in the code in case 
circumstances are found in which a smoother evolution 
is obtained. 



D. Implementation of the approximate 
conservation laws 



The theoretical basis for using approximate angular 
momentum and energy conservation to improve the accu- 
racy of the orbit computation, was discussed in Sec. lIII CI 
We now present details of how this is implemented for (1) 
the angular momentum in an equatorial orbit, (2) the 
angular momentum in a polar orbit, and (3) the energy. 
The code is written so that all of these approximate con- 
servation laws may be used, or not, simply by changing 
input parameter switches. 



1. Angular velocity in an equatorial orbit 

The angular momentum per unit mass 

h = qv p -pv q (4-11) 
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is approximately conserved. In terms of proper time r 
along the particle's trajectory, 

_ = u% + g ^W v? -p^£. (4.12) 

Now, Eq. (|4.1|l takes the form 

^v A = ~^((v p ) 2 + (v q ) 2 )+E A , A = (q,p), (4.13) 

where the Ea contain only small quantities. We also in- 
troduce the small quantity fj? — (g la —gJ2.^)v a , which rep- 
resents the difference between raising an index of the co- 
variant velocity by the full or background metric. Thus, 



A A V A 

v +72-- 



(4.14) 



Combining Eqs. (|4.12|) . (|4. lrif> and (|4.14|) . we obtain 



dh 



= fi q v p - [i v v q + qE p - pE q , 



(4.15) 



which is implemented in the code. We extract va from 
the evolved value of h. This is done by using the con- 
straint that the particle is on the equator, q 2 + p 2 = 1. 
Thus qv q + pv p = so that 



H/'" + ^J=0- (4.16) 
Combining Eqs. 14.11|l and l|4.16(l . we find 

v q = -ph-r 2 q(qfi q +pij, p ) 



v p = qh- r 2 p(qn q +pn p ), 



(4.17) 



which is implemented in the code. Furthermore, the 
particle is constrained to follow the equator exactly, 
and so the particle's position is corrected according to 
z A — > z A f e with 



fe 



1 



q z + p z 



(4.18) 



2. Angular velocity in a polar orbit 

In the case of polar motion, simplified here to the case 
p = 0, the equations analogous to Eqs. (|4.11|) . I|4.15|l and 
(|4~T7|) are 

Fv„ dh FA„ „ 2h 

h = -£ L > fr = ^r +qv ^ > Vq = T- (4 - 19) 



3. The energy 

The energy per unit mass v u is conserved along a 
geodesic in the Schwarzschild background. In this case 



VuS 



h 2 + r 2 +v 2 (r 2 -2Mr) 



2v r r 2 



(4.20) 



We take v u s, as defined above, to be an approximately 
conserved quantity. From Eq. 14. 1|) . 



d 



h 2 v 2 M 



Ex 



(4.21) 



Using 



dr r ( 2M\ 

— =v r = -v uS +(l- — )v r + f i r , (4.22) 



differentiation of Eq. (|4.20|l leads to 
dh 



d 

7? 



v uS = (2h—rv r + 2vlMrn r - 2h 2 v r \i T 
\ dr 

+E 1 {r 3 v 2 -r 3 - 2Mr 2 v 2 ~ h 2 r)) —J--. (4.23) 



There is an option in the code to evolve v u s by Eq. I|4.23|l . 
In this case, we extract v r from the value of v u s- This 
is done by rewriting Eq. (|4.20|) as a quadratic in v r . We 
find 



v u s ± 



'uS 



h 2 • 



(4.24) 



When the code is evolving v u s by Eq. (|4.23|) . at each 
time step it also evolves v r in the usual way. The ± in 
Eq. (|4.24() is chosen so that the result for v r is closest to 
the directly evolved value; further, if the square root in 
Eq. H4.24fl is less than some threshold, or imaginary, the 
directly evolved value of v r is not corrected. For a circular 
orbit of the Schwarzschild background, the square root is 
exactly zero, and therefore it is difficult to use this option 
when evolving a circular orbit. 



E. The metric variable W 

The only Bondi-Sachs metric variable that is non-zero 
in the background metric is W, and in order to improve 
numerical accuracy, the code treats W as the sum of 
the background analytic part (Wan) plus a correction 
(Wnum)- The values of Wan and its derivatives are found 
exactly, and finite differencing is applied only to the part 
Wnum- In effect, this also applies to the other metric 
variables, because their background analytic parts are 
zero. 



V. COMPUTATIONAL TESTS AND RESULTS 
A. Convergence 

1. Initial accelerations 

Firstly, we investigated the convergence of various ac- 
celerations on the initial null cone, and in so doing 
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tested the gravitational hypersurface equations, the grav- 
itational evolution equation and the particle evolution 
equations. The tests were made with the particle ini- 
tialized at r — 9 at the north pole with v r = v' p = 
and v q set to the value for a circular orbit. The parti- 
cle mass was m = 1CP 4 . The particle velocity was up- 
dated directly, without incorporating approximate con- 
servation laws. The overlap between north and south 
patches was minimal (q s = 1.0). The size of the poly- 
trope was = 5.0. The following quantities, all of 
which are rates of change, were determined on the ini- 
tial null cone: ||J u ||oo> h ;U v u<u , v TiT and v q , T . The 
quantities involving u— derivatives were found by evolv- 
ing the code for one time-step and then applying the 
formula Q. u = (Qi — Qo)/Au; the quantities involving 
r— derivatives are found directly by the code using data 
only on the initial null cone. 

The following grids were used: (a) coarse, n x = 41, 
n q = n p = 25; (b) medium n x = 81, n q = n p = 45; and 
(c) fine, n x = 161, n q — n p = 85. In the different grids, 
A x and A q = A p scale as 4:2:1. The (single) time-step 
was Au — 10~ 5 , which, for all grids, is much smaller than 
the spatial discretization, so that second order spatial ac- 
curacy is expected. Assuming that a quantity Q behaves 
as Q = a + bA n , it is straightforward to show that 



n = loe 



Qc Qm 
Qm Q f 



(5.1) 



where Q c , Q m and Qf refer to the computed values of Q 
using the coarse, medium and fine grids, respectively. 

Our results are stated in Table [I] it is clear that, on 
the initial null cone, the polytropic model is convergent 
with the order n in the range 1.59 to 2.28. 



2. Circular orbit 

Secondly, we performed a convergence test for a par- 
ticle in a circular orbit around a black hole. For the 
coarse and medium grids, the particle completed a whole 
orbit, but for the fine grid this was not possible. The 
particle was initialized at r — 9, q = and p = 1 
with v r — v p — and v q set to the value required for 
a circular equatorial orbit. The mass of the particle 
was m = 10 -6 and the size was i?* = 3 (the require- 
ment that the polytrope should be resolvable by all grids 
places a lower limit on i?»). We used the technique in 



TABLE I: Convergence of the Polytropic model 





Coarse 


Medium 


Fine 


n 




0.4346x10" 


2 


0.4441 xlO" 2 


0.4460xl0~ 2 


2.28 


h, u 


-0.7377x10" 


-2 


-0.3436xl0~ 2 


-0.2129xl0~ 2 


1.59 




0.2732 xl0~ 


4 


0.1273xl0~ 4 


0.0789xl0~ 4 


1.59 


Vr,r 


-0.5050x10" 


-3 


-0.5455xl0~ 3 


-0.5564xl0~ 3 


1.9 


Vq,T 


-1.8069x10" 


-3 


-0.8416xl0~ 3 


-0.5215xl0~ 3 


1.59 




FIG. 1: Convergence rate of E/m as a function of time for 
< u < 20. 



Sec. IIV Dl to model approximate conservation of angular 
momentum, but not of energy. The angular grid-patch 
overlap was q s — 1.2. The test results are for the time 
interval < u < 20 representing just under 1/8 of a 
whole orbit. The grids used were (a) coarse, n x = 61, 
n q — n p = 20 with du = 1.6666 x 10 -2 ; (b) medium, 
n x = 121, n q = n p = 35, du = 8.3333 x 10~ 3 ; and (c) 
fine, n x — 241, n q = n p = 65, with du = 4.1666 x 10~ 3 . 

The convergence rate of v u , between u = and u = 20 
as estimated from Eq. I|5.1|) . is shown in Fig. 2] the av- 
erage rate is n = 1.064, so that the effective convergence 
rate of the particle's energy is approximately first order. 
This is an artifact of the numerical scheme used in the 
evolution equation, which for a fixed value of the dissi- 
pation parameter is only first order in time pl|. 



B. Whole orbit with m/0 

The various input parameters were described in 
Sec. IV A 21 ( circular orbit) , medium grid. The computa- 
tion was run for 25,000 time steps until u = 208 and rep- 
resents more than one orbit (which is achieved at about 
u = 170); the computation took about 24 hours of wall- 
clock computer time. As discussed in Sec. IV A 21 the re- 
sults are within the convergence regime of the numerical 
method. The run was performed for illustrative purposes 
and is not physical because a polytrope with the param- 
eters used here would be tidally disrupted. 

The results of the computation are shown in Figs. [21 
to [SI in which the particle inspirals (Fig. |2J) , losing en- 
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FIG. 2: r-coordinate for a complete orbit. 



ergy (Fig. 0J and angular momentum (Fig. [3J • Figure 
shows the time development of the L\ norm of J. u , which 
provides a measure of the dynamic activity of the gravita- 
tional field. We also performed whole orbit computations 
with input parameters as above, but varying the particle 
mass by powers of 10 in the range 1CP 4 to 10 -9 . Wc 
found that the energy loss rate scales, as expected, as 
m 2 . 

In our numerical method, the measured inspiral was 
Ar = —0.016 after one complete orbit. The inspiral af- 
ter one orbit is Ar = —2.98 x 10 -6 according to the 
quadrupole formula [3(j, and Ar = —4.75 x 10 -6 ac- 
cording to a perturbative method [Til l37| . The rate 
of energy loss, i.e. the rate of change of E = mv u , is 
a better measure of the inspiral rate - because Ar af- 
ter one orbit includes a contribution due to the orbit 
becoming slightly elliptical. Averaged over a complete 
orbit, the measured rate for the numerical method is 
dE/du = 2.669 x 10~ 13 , whereas the quadrupole formula 
predicts dE/du = 1.10 x 10 -16 , and the perturbative 
method gives dE/du = 1.75 x 10 -16 . There is thus a dis- 
crepancy between the energy loss rates found here and 
by other methods. The cause of the discrepancy is not 
known, and may comprise a number of factors: (a) in 
order to resolve the particle properly, we are forced to 
make its size too large for it to be physical, and thus the 
model ignores internal tidal stresses that in this case are 
large (see Sec. lIII Aj) : (b) lack of resolution; and (c) other. 
The issue is discussed further in the Conclusion, Sec. lVIl 
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FIG. 3: Angular momentum per unit particle mass for a com- 
plete orbit. 
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FIG. 4: Energy per unit particle mass for a complete orbit. 



C. Capture of particle by the black hole 

The purpose of the test was to see how the code be- 
haves as the particle approaches the event horizon at 
r = 2, and the particle was initialized at r = 6, i.e. at 
the ISCO. The size of the particle was i?* = 2, and so, as 
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with the complete orbit computation, the model ignores 
significant tidal stresses and is not physical. In order to 
shorten the inspiral time the particle was given a small 
inward radial velocity (v r = —0.01), and the angular ve- 
locity was set to that for a circular orbit in the test par- 
ticle limit (v 9 = 0.0962, v p = 0). The test was performed 
with two different grids, enabling us to see at which stage 
numerical errors become significant. The grids used were 
n x = 121, n q = rip = 35, du = 8.3333 x 10~ 3 (medium 
grid); and n x = 81, n q — n p = 25, du — 1.25 x 10~ 2 
(coarser grid). 

In the coordinates being used, as r — > 2 the evolution 
variable v r — > oo. Thus, because of this coordinate ef- 
fect, we expect the code to crash at some value of r just 
greater than 2. The results of the computation are shown 
in Figs. H3to[51 In the medium grid computation, the 
particle inspirals until the code crashes at u = 182.5 with 
the particle at r = 2.00077 and \v r \ fts 5, 000. The parti- 
cle completed just over two complete revolutions, i.e. its 
angular position changed by just over Att radians during 
the evolution. The particle crossed r = 2.1 at u = 162.5, 
r = 2.01 at u = 172.5 and r = 2.001 at u = 181.7; 
thus demonstrating a freezing of radial position, as ex- 
pected due to the redshift inherent in the u-coordinate. 
Throughout the computation, the position of the parti- 
cle varies smoothly (Fig. 0). The particle loses energy 
(Fig. |SJ and angular momentum (Fig. [7J at a fairly con- 
stant rate, until about u — 150, r = 3.2. Further, the ac- 
tivity of the gravitational field (as measured by || J. u \\. 
Fig. starts to grow rapidly at this time. We have 
not analyzed the cause of this effect. The results of the 



FIG. 6: The orbit, in the (x,y) — (r cos (j>, r sin <j>) plane, 
traced by a particle, initially at r = 6, as it is captured by 
a black hole. Overlaid in the figure is the orbit traced by a 
particle initially at r — 9 (dotted line). The central circle 
indicates the location of the horizon (r = 2) . 
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FIG. 7: Angular momentum per unit particle mass for the 
capture of a particle by a black hole. 
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D. Other experiments and computational issues 

A run was performed using the same input parameters 
as in Sec. IV Bl except that m = 0. The orbit was found to 
be exactly circular, with no loss of energy or angular mo- 
mentum. Because the numerical method implements the 
approximate conservation of angular momentum, which 
in this case is exact, the angular momentum behaved as 
expected. However, conservation of energy was not en- 
forced. The only (minor) error was in the time for an 
orbit: the numerical method yielded Ait = 169.638 (for 
the grid used in this specific run), whereas the analytic 
value is 2nr^ = 169.646. In order to check that the ex- 
act conservation of energy was not just a consequence of 
the orbit being exactly circular, we also did a run for a 
non-circular orbit, starting at r = 9, q = , p = 1 and 
with initial velocity v r = 0, v q = 0.047, w p = 0. The 
particle was found to move in an orbit between r = 9 
and r — 11.79, and there was a very small change in v u : 
energy was conserved, in the sense of dv u /du = 0, to the 
order of one part in 10 14 . 

A simple test of the code is to determine if it produces 
results concerning caustic formation that are consistent 
with the estimate given in Eq. I|3.10fl . We performed a 
number of short runs (up to 100 iterations), each with 
a different value of particle mass to, and made a binary 
chop search to find the critical value of to above which 
the code crashes due to the onset of caustic formation 
(indicated by the metric variable (3 — > oo at T + ). The test 
was performed with initial velocity v % = 0, initial position 
r — 9 at the north pole, and polytropic radius i?„ = 2. 



-0.9428 



medium 

- - coarser 




FIG. 9: L\ norm of the radiation indicator, || J lU ||i, shown 
using a logarithmic scale in the vertical axis, for the capture 
of a particle by a black hole. 



The grid discretization was n x = 121, n q = n p = 35, 
with time step dt = 8.3333 x 10~ 3 . We found that the 
code behaved properly with to = 0.03 but crashed when 
to = 0.04, consistent with the critical value of m = 0.111 
indicated by Eq. I|3.10(l . 

Unfortunately, it is not possible, at this time, to 
present results on gravitational radiation output. The 
module used in the code for calculating the news [2l| 
was originally developed and tested under conditions in 
which the fields are well resolved at X + , which is not the 
case for particle applications. Improvements in the news 
computation have recently been investigated [38ll39l |. but 
it is not yet known whether the radiation from a particle 
source can be reliably computed. Results will be reported 
elsewhere after the necessary development and testing. 

The tests were performed on a Linux machine with a 
single processor running at 1.8GHz. The complete orbit 
run reported in Sec. IV Bl used a grid of 121 x 35 2 points 
(on each angular patch) and 25,000 time steps. The run 
time was about 24 hours. Of course, the run time scales 
with grid discretization as A . 



VI. CONCLUSION 



FIG. 8: Energy per unit particle mass for the capture of a 
particle by a black hole. 



In this paper we have described and implemented the 
PP method for evolving the full Einstein equations us- 
ing a characteristic evolution code. The method could 
be adapted and used within other numerical relativity 
frameworks. The PP method can be a useful tool in 
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modeling astrophysical situations involving a black hole 
and another much smaller object, in a regime in which 
a full hydrodynamic model is not necessary. We have 
demonstrated that the method works in the sense that 
it computes orbits that are convergent, and deviations 
from a Schwarzschild orbit scale as expected with par- 
ticle mass to. However, the computed inspiral rate was 
much larger than that predicted by other methods. A 
feature of the code is that it avoids the computational 
expense of modeling the material object by means of rel- 
ativistic hydrodynamics. 

The PP method has the potential to supply accurate 
orbits, including inspirals towards the ISCO and plunges 
to the black hole, as well as the associated gravitational 
radiation output. Such results would be useful both di- 
rectly, and also indirectly in providing error bars against 
which other methods could be tested. In the tests de- 
scribed here, the need to resolve the particle has pre- 
vented us from making the particle radius small enough 
so that its local rigidity is justified. In order to achieve 
a proper physical basis for the model, we envisage the 
following future work 

1. It is necessary to investigate the effect of the poly- 
trope radius (i?*) on particle motion. On physi- 
cal grounds, one would expect that in some regime 
the particle motion should be independent of i?*. 
(Of course, taking the limit i?* — s- is equivalent 
to changing from the polytropic to a d— function 
model, and this link provided motivation for the 
investigation of that model) . 

2. Once parameters for the model that are physically 
realistic can be attained, it will be necessary to in- 
vestigate whether the PP method computes reliable 
energy loss rates. 

3. The gravitational radiation output (Bondi news 
function) is required, both to supply a waveform 
and to check the energy balance (the rate of loss of 

I 



orbital energy mv u should be of the same magni- 
tude as the radiation power). 

4. Once the above issues have been resolved, it will 
be necessary to validate the PP method by obtain- 
ing results that agree, in some regime, with results 
obtained by another method. 
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APPENDIX A: THE GEODESIC EQUATION 

We have used Maple to compute the form of Eq. I|4.1|) 
for the metric (|2.1|l . The angular part of Vi, va is rep- 
resented by the spin weighted quantity V ang = vaQ A - 
Further, for ease of application to the approximate con- 
servation method fSecs. ITTlCI and HVDfl . the formulas 
are presented with the zeroth order quantities (in each 
case, the first line) shown separately from the perturba- 
tive (Ei) quantities. 



dvr_ = 2V ang V ang -v?(rW, r -W)r 
dr 2r 3 

- rV ang V ang JJ r - ^ r e- 2 ^v r r 3 UV ang K - <r e-^v r r 3 UV ang K - 2JV 2 ng K 

+ 4e- 2/ V/3 r (r + W)r 2 K - 80 >r e- 20 v r r 3 v u K + r.J, r V 2 ng K + 2e~ 20 v r r 3 U, r V ang K 

+ 2e-^v r r 3 U r V ang K + rJ r V 2 ng K^j^- E ; (Al) 

dV ang _ (g + ip)VgngVang 

dr r 2 
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+ ^ - % + ip)V„ nfl - 1) + 4(9^) JJK„ s Kn S - 2(gj)XJK„ ff K„ s - 2(SJ)tfJV cms V cmg 

- 4e- 2 Pv r (d[3)r 2 UV ang - 4e- 2 V(3/3V 2 £7K„ s - 2e~^ v 2 r {^W)r + ^ 2Si v 2 r (r + W)(8/3)r 

-- 8e- 2 M3/?>V + (3J)J 2 K 2 „ S + (dJ)j 2 V 2 ng - 2(<5K)JKV 2 ng - 2(dK)JKV 2 ng + 2(dK)V ang V ang 

+ (SJ)V a 2 ng JJ + mV 2 ng JJ + 4ipJV 2 ng + (dJ)V 2 ng + (dJ)V 2 ng + 4qJV 2 ng + 2r 2 e- 2 f } v r (W)V ang 

+ 2r 2 e- 2t3 v r (W)V ang + Ar 2 e- 2fi v r qUV ang + Air 2 v rP UV ang ) (A2) 
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